function [X1, Y1, X2, Y2] = suh_radcliffe(X0, Y0, X0_1, Y0_1, x1, y1, x2, y2, x3, y3, th12, th13)
    % SUH_RADCLIFFE calculates the coordinates of the movable joints (X1, Y1) and (X2, Y2)
    % in a four-bar linkage mechanism using the displacement matrix method.
    % 
    % References:
    % Suh, C. H. and Radcliffe, C. W., "Synthesis of Plane Linkages With Use of the Displacement Matrix," 
    % Journal of Engineering for Industry, vol. 89, no. 2, pp. 206-214, 1967. DOI: 10.1115/1.3610029
    %
    % Input:
    %   - X0, Y0: Coordinates of the first fixed hinge
    %   - X0_1, Y0_1: Coordinates of the second fixed hinge
    %   - x1, y1: Coordinates of the homologous point in the first configuration
    %   - x2, y2: Coordinates of the homologous point in the second configuration
    %   - x3, y3: Coordinates of the homologous point in the third configuration
    %   - th12: Rotation angle between configuration 1 and 2 (radians)
    %   - th13: Rotation angle between configuration 1 and 3 (radians)
    %
    % Output:
    %   - X1, Y1: Coordinates of the first movable joint
    %   - X2, Y2: Coordinates of the second movable joint
    %
    % Author: Lorenzo De Sanctis

    %% Constants
    R2 = x2 - x1 * cos(th12) + y1 * sin(th12);
    S2 = y2 - x1 * sin(th12) - y1 * cos(th12);
    R3 = x3 - x1 * cos(th13) + y1 * sin(th13);
    S3 = y3 - x1 * sin(th13) - y1 * cos(th13);

    %% First movable joint
    % Position 1-2
    A12 = R2 * cos(th12) + S2 * sin(th12) - X0 * cos(th12) - Y0 * sin(th12) + X0;
    B12 = S2 * cos(th12) - R2 * sin(th12) + X0 * sin(th12) - Y0 * cos(th12) + Y0;
    C12 = R2 * X0 + S2 * Y0 - 0.5 * (R2^2 + S2^2);

    % Position 1-3
    A13 = R3 * cos(th13) + S3 * sin(th13) - X0 * cos(th13) - Y0 * sin(th13) + X0;
    B13 = S3 * cos(th13) - R3 * sin(th13) + X0 * sin(th13) - Y0 * cos(th13) + Y0;
    C13 = R3 * X0 + S3 * Y0 - 0.5 * (R3^2 + S3^2);

    %% Second movable joint
    % Position 1-2
    A12_1 = R2 * cos(th12) + S2 * sin(th12) - X0_1 * cos(th12) - Y0_1 * sin(th12) + X0_1;
    B12_1 = S2 * cos(th12) - R2 * sin(th12) + X0_1 * sin(th12) - Y0_1 * cos(th12) + Y0_1;
    C12_1 = R2 * X0_1 + S2 * Y0_1 - 0.5 * (R2^2 + S2^2);

    % Position 1-3
    A13_1 = R3 * cos(th13) + S3 * sin(th13) - X0_1 * cos(th13) - Y0_1 * sin(th13) + X0_1;
    B13_1 = S3 * cos(th13) - R3 * sin(th13) + X0_1 * sin(th13) - Y0_1 * cos(th13) + Y0_1;
    C13_1 = R3 * X0_1 + S3 * Y0_1 - 0.5 * (R3^2 + S3^2);

    %% Linear system
    A = [A12, B12, 0, 0; 
        A13, B13, 0, 0; 
        0, 0, A12_1, B12_1;
        0, 0, A13_1, B13_1];
    b = [C12; C13; C12_1; C13_1];
    x = A \ b;

    % Extract coordinates
    X1 = x(1);
    Y1 = x(2);
    X2 = x(3);
    Y2 = x(4);
end